#Pkg.add("WAV")
#Pkg.add("Plots")
using WAV,Plots
plotly()
Im Folgenden arbeiten wir mit zwei Klangdateien.
(s1,fs1) = wavread("../speech1.wav");
(s2,fs2) = wavread("../speech2.wav");
println("Abtastrate für speech1.wav: ",fs1," Hz")
println("Abtastrate für speech2.wav: ",fs2," Hz")
If s is a signal with sampling frequency fs, then the i-th measurement in s takes place at time (i-1)*period_length where period_length = 1/fs.
(Note: Indexing of vectors starts at 1 in both Julia and Matlab)
v_sample_time1 = (0:(length(s1)-1))/fs1;
v_sample_time2 = (0:(length(s2)-1))/fs2;
plot(v_sample_time1,s1;title="speech1.wav",xlab="Zeit [s]")
plot(v_sample_time2,s2;title="speech2.wav",xlab="Zeit [s]")
TODO
Gegeben sind
v_signalsampling_rateframe_lengthframe_shift (uff, das geht sicher eloquenter)Die erste Aufgabe besteht darin die Abmessungen für unsere Zeitfenster, die in Sekunden angegeben sind, als eine Anzahl von Samples auszudrücken. Dadurch ergeben sich folgende Definitionen:
samples_per_frame = Int(floor(frame_length * sampling_rate))samples_per_shift = Int(floor(frame_shift * sampling_rate))Als nächstes stellt sich die Frage wie viele Fenster in unser Signal passen.
Wir nehmen sinnvollerweise an, dass in unser Signal mindestens ein Fenster (genauer gesagt, das erste Fenster) passt.
Jedes darauf folgende Fenster konsumiert von den restlichen length(v_signal) - samples_per_frame Samples genau samples_per_shift.
Dementsprechend ergibt sich die Formel
shift_count = 1 + Int(floor((length(v_signal) - samples_per_frame)/samples_per_shift))Diese Daten reichen nun aus, um die Ausgabematrix zu befüllen.
Für die Mittelpunkte der Zeitfenster ergibt sich durch folgende Überlegung: Der Mittelpunkt des ersten Fensters befindet sich bei frame_length/2. Jeder darauffolgende Mittelpunkt ergibt sich aus einer Verschiebung um frame_shift.
v_time_frame = frame_length/2 + frame_shift * (0:(shift_count-1))my_windowing = function(v_signal, sampling_rate, frame_length, frame_shift)
samples_per_frame = Int(floor(frame_length * sampling_rate))
samples_per_shift = Int(floor(frame_shift * sampling_rate))
shift_count = Int(floor((length(v_signal) - (samples_per_frame - samples_per_shift))/samples_per_shift))
m_frames = zeros(samples_per_shift, shift_count)
for j in 1:shift_count
for i in 1:samples_per_shift
m_frames[i,j] = v_signal[samples_per_shift*(j-1)+i]
end
end
v_time_frame = frame_length/2 + frame_shift * (0:(shift_count-1))
(m_frames, v_time_frame)
end
(frames1,frame_times1) = my_windowing(s1,fs1,32e-3,16e-3);
(frames2,frame_times2) = my_windowing(s2,fs2,32e-3,16e-3);
autocor = x -> conv(x,reverse(x))/length(x)
truncated version of autocor
autocor_t = x -> autocor(x)[1:length(x)]
ff_estimator = function(s,fs)
(frames, frame_times) = my_windowing(s,fs,32e-3,16e-3)
(samples_per_frame, shift_count) = size(frames)
ac = zeros(frames) # auto correlations
freq = zeros(frame_times) # ff estimate per frame
# 80Hz:400Hz frequency window
lb = Int(floor(fs/400)) # lower bound
ub = Int(floor(fs/80)) # upper bound
for j in 1:shift_count
ac[:,j] = autocor_t(frames[:,j])
period = lb - 1 + indmax(ac[lb:ub,j]) # period estimate (measured in number of samples)
freq[j] = fs/period
end
return (frame_times, freq)
end
Signal strengths s1 and fundamental frequency estimates ff1 have different units (and magnitudes) which is why we have to rescale in order to see the details of both data in one plot.
(tf1,ff1) = ff_estimator(s1,fs1);
plot([v_sample_time1,tf1],[s1*2000,ff1])
(tf2,ff2) = ff_estimator(s2,fs2);
plot([v_sample_time2,tf2],[s2*2000,ff2])
TODO